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Abstract 

We study numerically statistical distributions of sums of chaotic orbit coordi- 
nates, viewed as independent random variables, in weakly chaotic regimes of 
three multi-dimensional Hamiltonian systems: Two Fermi-Pasta-Ulam (FPU- 
/3) oscillator chains with different boundary conditions and numbers of particles 
and a microplasma of identical ions confined in a Penning trap and repelled by 
mutual Coulomb interactions. For the FPU systems we show that, when chaos is 
limited within "small size" phase space regions, statistical distributions of sums 
of chaotic variables are well approximated for surprisingly long times (typically 
up to t ~ 10 6 ) by a g-Gaussian (1 < q < 3) distribution and tend to a Gaussian 
(g=l) for longer times, as the orbits eventually enter into "large size" chaotic 
domains. However, in agreement with other studies, we find in certain cases that 
the g-Gaussian is not the only possible distribution that can fit the data, as our 
sums may be better approximated by a different so-called "crossover" function 
attributed to finite-size effects. In the case of the microplasma Hamiltonian, 
we make use of these g-Gaussian distributions to identify two energy regimes 
of "weak chaos" -one where the system melts and one where it transforms from 
liquid to a gas state-by observing where the g-index of the distribution increases 
significantly above the q = 1 value of strong chaos. 

Keywords: Nonextensive Statistical Mechanics, g-Gaussian distributions, 
"edge of chaos", weak chaos, quasi-stationary states, multi-dimensional 
Hamiltonian systems 
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1. Introduction 



Probability density functions (pdfs) of chaotic trajectories of dynamical sys- 
tems have been studied for many decades and by many authors, aiming to un- 
derstand the transition from deterministic to stochastic dynamics [li-lllf. The 
fundamental question in this context concerns the existence of an appropri- 
ate invariant probability density (or ergodic measure), characterizing chaotic 
motion in phase space regions where solutions generically exhibit exponential 
divergence of nearby trajectories. If such an invariant measure can be estab- 
lished for almost all initial conditions (i.e. except for a set of measure zero), one 
has a firm basis for studying the system from a Statistical Mechanics point of 
view. 

If, additionally, this invariant measure turns out to be a continuous and 
sufficiently smooth function of the phase space coordinates, one can invoke the 
Boltzmann-Gibbs (BG) microcanonical ensemble and attempt to evaluate all 
relevant quantities of equilibrium Statistical Mechanics, like partition function, 
free energy, entropy, etc. of the system. On the other hand, if the measure is 
absolutely continuous (as e.g. in the case of the so-called Axiom A dynamical 
systems), one might still be able to use the formalism of modern er godi c theory 
and SRB measures to study the statistical properties of the model [lCf . 

In all these cases, viewing the values of one (or a linear combination) of com- 
ponents of a chaotic solution at discrete times t n , n = 1, . . . , Af as realizations of 
Af independent and identically distributed random variables X n and calculat- 
ing the distribution of their sums, one expects to find a Gaussian distribution, 
whose mean and variance are those of the X„'s. This is indeed what happens 
for many chaotic dynamical systems studied to date which are ergodic, i.e. al- 
most all their orbits (except for a set of measure zero) pass arbitrarily close to 
any point of the constant energy manifold, after sufficiently long enough times. 
What is also true is that, in these cases, at least one Lyapunov exponent is 
positive, stable periodic orbits are absent and the constant energy manifold is 
covered uniformly by chaotic orbits, for all but a (Lebesgue) measure zero set 
of initial conditions. 

But then, what about "small size" chaotic regions of Hamiltonian systems, 
at energies where the maximal Lyapunov exponent is positive (but still rather 
"small") and stable periodic orbits still exist, whose islands of invariant tori and 
sets of cantori surrounding them occupy a positive measure subset of the energy 
manifold? In such regimes of so-called "weak chaos" , a great number of orbits 
stick for long times to the boundaries of these islands and chaotic trajectories 
diffuse slowly through multiply connected regions in a highly non-uniform way 
12 Such examples occur in many physically realistic systems studied in 



the current literature (see e.g. [15H18J). 

In this paper, we focus on such "weakly chaotic" regimes and demonstrat e by 
means of numerical experiments in the spirit of the Central Limit Theorem j!9| 
that pdfs of sums of orbit components do not rapidly converge to a Gaussian, but 
are well approximated, for long integration times, by the so-called q-Gaussian 
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distribution [2(| 



P(s) — aexp g (—j3s ) = a 



1 - (1 - q)Ps 2 



(1) 



where the index satisfies 1 < q < 3, /3 is an arbitrary parameter and a is 
a normalization constant. At longer times, of course, chaotic orbits eventually 
seep out from smaller regions to larger chaotic seas, where obstruction by islands 
and cantori is less dominant and the dynamics is more uniformly ergodic. This 
transition is signaled by the g-index of the distribution ([1]) decreasing towards 
q = 1, which represents the limit at which it becomes a Gaussian. 

Thus, in our models, ^-Gaussian distributions represent quasi- stationary 
states (QSS) that are often very long-lived, especially inside "thin" chaotic layers 
near periodic orbits that have just turned unstable. This suggests that it might 
be useful to study these pdfs (as well as their associated q values) for sufficiently 
long times and try to derive useful information about the corresponding QSS di- 
rectly from the chaotic orbits of the system, at least for time intervals accessible 
by numerical integration. QSS have already been studied in coupled standard 



maps and the so-called HMF model in [2l|, |22| . but not from the viewpoint 
of sum distributions. These authors calculated the growth of the variance of 
the total "angular momentum" of their models and obtained interesting results, 
which we shall discuss in the Conclusions section of the present paper. 

Care must be taken, however, with regard to the kind of functions that 
are used to approximate these pdfs. While it is true that g-Gaussians offer a 
popular and quite accurate representation of QSS sum distributions, they are 
not the only possible choice. As we have found in certain cases, other so-called 



"crossover" functions [20, |23|, |24( (see e.g. Eq. (TT4|) below) can be used to 
describe the same data with better accuracy. Thus, as pointed out already by 
several authors there may be other alternatives to q-Gaussians, that 



can approximate better sum distributions of weakly chaotic QSS such as the 
ones studied here. 

In what follows, after some preliminaries about our method in Section [21 we 
begin in Section [3] by studying the solutions of a Fermi-Pasta- Ulam (FPU)-/? 
1-dimensional chain of N = 128 particles, under periodic boundary conditions, 
near one of its simple periodic orbits called the 7r-mode, at energies where it 
has just turned unstable. This is one of the nonlinear normal modes (NNMs) of 
the chain and represents a continuation of the linear mode that consists of all 
neighboring pairs of particles oscillating with the same amplitude and a 7r-phase 



difference with respect to each other [281-13 



Following the approach of [34j , we verify that the pdfs of chaotic orbits start- 
ing very close to the unstable 7r-mode are well approximated by a q-Gaussian 
over time intervals of the order of t « 10 6 . These trajectories, however, represent 
a QSS: Their distributions for longer times deviate from a q-Gaussian and show 
a tendency to converge to a Gaussian as q — > 1. Remarkably enough, this result 
holds true both for the distribution of the single values of the chaotic variable 
(as shown recently in [34|). as well as the distribution of its sums relevant to 
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our study. 

In seeking to relate our statistical results to the actual dynamics of the chain, 
we adopt a very interesting approach described by Cretegny et al. in [35| and 
compute their Co(t) function for initial conditions close to the 7r-mode at low 
energies, just above the destabilization threshold of the mode (always under 



periodic boundary conditions). We find, just as in [35[, that Co(t) displays a 
maximum as a chaotic breather appears on the chain, until t ss 10 s , when the 
breather collapses and energy equipartition occurs. Remarkably, this is about 
the time when our q-Gaussian approximations converge to Gaussians. On the 
other hand, for an FPU chain with ./V = 129 particles under fixed boundary 
conditions the Co (t) function does not exhibit a clear maximum and no chaotic 
breather is observed, while a q-Gaussian approximation is still valid for long 
times, eventually turning into a Gaussian when Co(t) settles down to its limiting 
value and energy is equally distributed among all degrees of freedom. 

We then turn, in Section [51 to a similar investigation regarding an FPU-/3 
model with only N — 5 particles and fixed boundary conditions. We first choose 
a NNM called the SPOl mode where, between every two particles oscillating 
with the same amplitude and opposite phase, there is one that is stationary for 
all t. Starting again with a total energy just above the (first) destabilization 
threshold and choosing initial conditions in the vicinity of the NNM, we discover 
that the sum distributions of chaotic orbits have again a (/-Gaussian shape, with 
q significantly higher than 1. Integrating the equations of motion for times of 
the order of t = 10 6 , we observe that the orbits form a thin "figure eight" in a 
2-dimensional (91, pi) projection of the energy manifold. In fact, the closer one 
starts to the NNM, the longer the persistence of the (/-Gaussian QSS and the 
thinner the width of the "figure eight" region. 

However, as soon as we start deviating significantly from the unstable pe- 
riodic point, the transitory nature of the QSS becomes apparent: Already for 
time intervals of the order of t = 10 5 or t = 10 6 , sum distributions begin to show 
a clear tendency to approach a Gaussian. Furthermore, the "figure eight" -like 
region expands to larger domains of phase space, where the Lyapunov exponents 
increase by orders of magnitude and the motion is more "uniformly" chaotic. 
Analogous results in this direction have also been obtained very recently on 
low-dimensional conservative maps in 



We then carry out a similar study near another NNM (i.e. the so-called 
SP02 mode), where every two oppositely moving particles one is stationary for 
all times. In this case we find QSS that are much longer-lived, possessing sum 
pdfs of the q-Gaussian type that do not tend to a Gaussian even for times as 
long as t = 10 10 ! Here, chaotic orbits trace out thin regions of nearly vanishing 
Lyapunov exponents close to banana-shaped tori representing a true "edge of 
chaos" regime. The interesting observation here is that the pdf we obtain in 
the long time limit converges to a smooth function that is well approximated by 
the "crossover" formula (fl"4")l proposed to hold in a regime between q-Gaussians 
and Gaussians, where finite size (and time) effects yield distributions with lower 
tails [13, [H, [24 1 . 

Finally, we turn in Section [6] to a very different Hamiltonian system that 
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describes N identical ions in a microplasma, mutually repelled by Coulomb 
interactions and confined radially by a magnetic field and axially by an elec- 
trostatic field (Penning trap). Here the interactions are long range, unlike the 
FPU models where only nearest neighbors are involved. In a recent publication 



371 ] . we have shown that there is a critical energy range, 2 < E < 4, where 



N = 5 ions undergo a "melting transition" [38], as they exit the domain of 
their local potential minima and begin to wander chaotically throughout the 
available phase space. Interestingly enough, when one uses g-Gaussian distri- 
butions to approximate the statistics of orbits in this region, one discovers that, 
for 2 < E < 3, the q values steadily increase above I and attain a peak value 
of q ps 1.8 (near E pa 2.8) before returning to q = 1 at higher energies of much 
larger positive Lyapunov exponents and uniformly spread chaos. Performing 
then a similar analysis for E > 3, we identify a second much broader energy 
range (i.e. E g (30,200)) over which the system undergoes a further transi- 
tion, passing from a liquid to a gas state, where q values attain another peak 
well above unity after which they start gradually to decrease and then fluctuate 
around the q = 1 value of uniform chaos. 

Regarding the computation of the Lyapunov exponents, we have applied 



throughout this work the well-known method of Benettin et. al |3jj|40j], together 
with a strategy proposed in 41| to follow the time evolution of the deviation 
vectors. Finally, we conclude with a summary and discussion of our results in 
Section [7] 



2. Statistical distributions of chaotic QSS and their computation 

The problems we investigate in this paper are described by an autonomous 
N degree of freedom Hamiltonian function of the form 

H = H(q(t), p(t)) = H( qi (t), . . .,q N (t), Pl (t), . . . ,p N (t)) = E (2) 

where (qk(t),Pk(t)) are the positions and momenta respectively representing 
the solution in phase space at time t. As is well-known, these solutions can be 
periodic, quasi-periodic or chaotic depending on the initial condition (q(0),p(0)) 
and the value of the total energy E. What we wish to study here are the 
statistical properties of such systems in reg imes of "weakly" chaotic motion, 
where the Lyapunov exponents [l(J Hil, Ho, 42 1 are positive but very small. 



Such situations often arise when one considers solutions which slowly diffuse 
into thin chaotic layers and wander through a complicated network of higher 
order resonances, often sticking for very long times to the boundaries of islands 
constituting the so-called "edge of chaos" regime [20j . 

There are several interesting questions to ask here: How long do these 
"weakly" chaotic states last? Assuming they are quasi-stationary, can we de- 
scribe them statistically by simply observing some of their chaotic orbits? What 
type of distributions characterize these QSS and how could one connect them to 
the actual dynamics of the solutions of the corresponding Hamiltonian system? 
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Does the dimensionality of the Hamiltonian, i.e. the number N of degrees of 
freedom, play a role in these considerations? 

The approach we shall follow is in the spirit of the well-known Central Limit 
Theorem [19f . In particular, we shall use the solutions of Hamilton's equations 
of motion 

dq k _ dH dp k dH _ 

dt dp k ' dt dq k ' 

to construct distributions of suitably rescaled sums of M values of a generic 
observable rji = rj(ti) (i = 1, . . . , M) which depends linearly on the components 
of the solution. If these are viewed as independent and identically distributed 
random variables (in the limit M — > oo), we may evaluate their sum 

M 

^=Y,^ ( 4 ) 

for j = 1, . . . , ATj c different initial conditions. Thus, we can study the statistics of 

the variables S$ , centered about their mean value (Sj9) — jj— 2i=i Vi^ 

ic J 

and rescaled by their standard deviation cxm 




(5) 



where 



M N- 

1C j=l 



^ E( 5 m - ^)) 2 = (^ )2 ) - (^) 2 - (6) 



Plotting the normalized histogram of the probabilities P{s < f)) as a function of 

(7) 

, we then compare our pdfs with a g-Gaussian of the form 



P(s$) = aexp 9 (-/34 )2 ) = a 



1 - (1 - q)/3.° m 



M 



(7) 



where q is the so-called entropic index, j3 is a free parameter and a a normaliza- 
tion constant [20( ■ Function ([7]) is a generalization of the well-known Gaussian 
distribution, since in the limit q — > 1 we have exp q (—(3x 2 ) — > cxp(—f3x 2 ). 

Moreover, it can be shown that the g-Gaussian distribution ([7]) is normalized 
when 

P = aV* '^7) (8) 

where T is the Euler T function. Clearly, Eq. d5J) shows that the allowed values 
of q are 1 < q < 3. 



6 



The index q appearing in Eq. ([7J is connected with the Tsallis entropy [2(| 



-, _ v^W q w 

S = k ^= lPl with V Pl = 1 (9) 

where i = 1, . . . , W counts the microstates of the system, each occurring with 
a probability pi and k is the so-called Boltzmann universal constant. Just as 
the Gaussian distribution represents an extremal of the BG entropy Sgg = 
Si = k Yli=x Pi m Pii so i s the g-Gaussian ([7]) derived by optimizing the Tsallis 
entropy ([9]) under appropriate constraints. 

Systems characterized by the Tsallis entropy are said to lie at the "edge of 
chaos" and are significantly different than BG systems, in the sense that their 
entropy is nonadditive and gene rally nonextensive 2(| 24\. In fact, a ^-Central 



Limit Theorem was proved [43| for g-Gaussian distributions ([7]). However, the 
role of this theorem in establishing the theoretical foundation of such functions 
as limiting distributions of systems of strongly correlated variables has been 
recently challenged in [27j |. 

Let us now describe the numerical approach we are going to follow to calcu- 
late the above pdfs. First, we specify an observable denoted by r](t) as one (or 
a linear combination) of the components of the position vector q(i) of a chaotic 
solution of Eq. ([3|), located initially at (q(0),p(0)). 

Assuming that the orbit visits all parts of a QSS during the integration 
interval < t < if, we divide it into N^ c equally spaced, consecutive time 
windows, which are long enough to contain a significant part of the orbit. Next, 
we subdivide each such window into a number M of equally spaced subintervals 
and calculate the sum S$ of the values of the observable r](t) at the right edges 
of these subintervals (see Eq. (|4])). 

In this way, we treat the point at the beginning of every time window as a 
new initial condition and repeat this process JVj c times to obtain as many sums 
as required for reliable statistics. Consequently, at the end of the integration, 
we compute the average and standard deviation of the sums dU, evaluate the 
N- K rescaled quantities s^) and plot the histogram P(s$) of their distribution. 

As we shall see in the next sections, in regions of "weak chaos" these dis- 
tributions are well-fitted by a q-Gaussian of the form J7j for fairly long time 
intervals. However, for longer times, as the orbits begin to diffuse to wider 
domains of "strong chaos" , the well-known form of a Gaussian pdf is recovered. 



3. FPU 7r-mode under periodic boundary conditions 

Let us start by considering the famous I-dimensional lattice of N particles 
with nearest neighbor interactions governed by the FPU-/? Hamiltonian [3, |45| 



N N 

tf(q,p)4E^ 2 + E 



lit 



E 



(10) 
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where qj is the displacement of the jth particle from its equilibrium position, 
Pj is the corresponding momentum, /3 is a real positive constant and E is the 
constant energy of the system. 

We first turn our attention to the well-known 7r-mode of Eq . 
has been studied in detail in several publications 28, 
defined by 

q.j{t) = -q j+1 (t) = q(t), j = 1, . . . , N (11) 
and exists only when periodic boundary conditions are imposed 

q N+k (t) = q k {t), k= 1,...,N; Vt (12) 

with N even (for a comprehensive review of NNMs and their connection with 
discrete symmetries of the Hamiltonian, see [HI). 

Our aim is to investigate here a chaotic region near this orbit at energies 
where it has just become unstable. Thus, we choose as our observable the 
quantity 

V(t)=q f (t) + q f _ 1 (t) (13) 

which satisfies T](t) — exactly at the 7r-mode and remains close to zero at 
energies E where the 7r-mode is still stable. However, r](t) will eventually deviate 
from zero at energies above the first destabilization threshold, i.e. E > E\. To 
compare our results with those published recently in [34 ] . we restrict ourselves 
to the case of N = 128 and j3 = 1 for which E\ w 0.0257 [12] and take as 
our total energy E = 0.768 (i.e. e = E/N = 0.006), at which our 7r-mode is 
certainly unstable. 

Regarding the accuracy of our integration, it is determined at each time step 
by requiring \H (q(t) , p(t)) — £j n | < 10~ 5 , where (q(t),p(t)) is the numerical 
solution obtained using the Burlisch-Stoer (non-symplectic) integrator with an 
adaptive step size control [46[ and i? m is the initial energy at time t = 0. We 
note that we have compared this method to Yoshida's 4th order symplectic 
integrator [ItJ with time step equal to 0.05 (yielding a relative energy error of 
order 10~ 6 - 10~ 7 ) and have found nearly identical results in every case. 

As we see in Fig. [TJ when we increase the total integration time if, the pdfs 
(red curves) we obtain approach closer and closer to a Gaussian with q tending 
to 1. Moreover, this seems to be independent of the values of iVj c and/or M, 
at least up to the final integration time tf = 10 8 that we have been able to 
check. For example, when we vary in Fig. QJb) the parameters iVj c and M, we 
obtain q-Gaussians of very similar shape with q between 1.51 and 1.67. It is 
important to note, however, that the same data may be better fitted by other 
similar looking functions: For example, we have found in the case of Fig. [TJa) 
that the numerical distribution (red curve) is more accurately approximated by 



a function presented in [20J, |23j, \2A \ 



P{s ( m) = — , a 1: a q > and q > 1 (14) 



{l-^ + ^exp[( g -l) ai 4f]} ; 
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Figure 1: Plot in linear-log scale of numerical (red curve), g-Gaussian (green curve) and 
Gaussian (blue curve) distributions for the FPU 7r-mode with periodic boundary conditions 
for N = 128, /3 = 1 and e = 0.006 > e u ~ 0.0002. Panel (a) corresponds to final integration 
time tf = 10 5 using JVj c = 10 4 time windows and M = 10 terms in the computation of the 
sums. In this case the numerical fitting with a q-Gaussian gives q » 1.818 with x 2 ~ 0.0007. 
Panel (b) corresponds to = 10 6 , iV; c = 10 4 and M = 100. Here, the numerical fitting gives 
q RJ 1.531 with x 2 ~ 0.0004. Panel (c) corresponds to t f = 10 s , N ic = 10 5 and M = 1000. 
It is evident that the numerical distribution (red curve) has almost converged to a Gaussian 
(blue curve). Panel (d) shows the same red pdf as in panel (a) together with the P function 
of Eq. Q3} for oi « 0.009, a q ss 2.849 and q « 2.179 with \ 2 ~ 0.000 08 (green curve). 



where a x w 0.009, a q w 2.849 and q w 2.179 with x 2 ~ 0.000 08, in contrast 
to the x 2 ~ 0.0007 obtained by fitting the same distribution by a g-Gaussian 
with q w 1.818 (see Fig. HJa)). Eq. (|14p represents the crossover between q- 
Gaussians and Gaussians and takes into account finite size (and time) effects 
reflected in the lowering of the tails of our distributions. 

Thus, we conclude that the chaotic state we have obtained is a QSS. Indeed, 
when we increase the final integration time beyond £f ~ 4 x 10 7 , we observe that 
the Lyapunov exponents monotonically increase and attain bigger values than 
those computed up to tt ~ 4 x 10 7 (see panel (b) of Fig. [3]). This may signify 
that the trajectories drift away from the neighborhood of the 7r-mode and enter 
a larger chaotic subspace of the energy manifold. Indeed, calculating a function 
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D(t) representing the euclidean distance of our orbit from the exact 7r-mode, 
we observed a sharp rise in the D(t) values at t f=s 7.5 x 10 8 . 





-30 -20 -10 10 20 30 




Figure 2: Plot in linear-log scale of numerical (red curve), (/-Gaussian (green curve) and 
Gaussian (blue curve) distributions for the FPU 7r-mode with periodic boundary conditions 
for N = 128, P = 1 and e = 0.006 > e u ~ 0.0002 following [3^|. £ is defined as in [3 as well. 
Panel (a) corresponds to a final integration time of tf = 10 6 . In this case the numerical fitting 
with a (/-Gaussian gives q m 1.495 with x 2 ~ 2.63 X 10 -5 . Panel (b) corresponds to t { = 10 7 . 
Here, an analogous fitting gives q Ri 1.374 with \ 2 ~ 3.86 X 10~ 5 . Panel (c) corresponds to 
if = 10 8 getting q w 1.279 with x 2 ~ 9-86 x 10 -6 . 



Similar results concerning the tendency of the distributions to approach a 
Gaussian for > 10 6 are also obtained using the methodology proposed recently 



34j . In this paper, however, the authors use a different way of calculating 



the numerical distributions of the rji observables: Instead of using sums, they 
divide their total time interval [0, tf] into short segments of 100 integration steps 
each and study the pdfs of their observables evaluated at the right endpoints 
of these segments. Moreover, since the longest integration time they considered 
was if = 10 6 , they only observed QSS approximated by g-Gaussians and did 
not detect their transitory character, which becomes evident for longer times. 
Indeed, when one takes if = 10 7 or 10 8 , one finds that distributions computed 
by the strategy of [34[ , show a tendency to approach a Gaussian by finding that 
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q — > 1 as we can see in Fig. [5] of the present paper. 



Chaotic breathers and the FPU 7r-mode under Periodic Boundary 
Conditions 



In [35J , the authors studied the time evolution of an FPU-/3 chain towards 
equipartition by considering initial conditions close to the 7r-mode under periodic 
boundary conditions. They observed that at energies well above the threshold 
where this mode becomes modulationally unstable, a remarkable localization 
phenomenon occurs. In particular, they report the spontaneous appearance 
of excitations that strongly resemble discrete breather solutions, but have a 
finite lifetime, while their dynamics is chaotic as indicated by their positive 
Lyapunov exponents. Moreover, their numerical results suggest that the lifetime 
of these chaotic breathers is related to the time necessary for the system to reach 
equipartition. 

The appearance of chaotic breathers produced as time evolves can be mon- 
itored by evaluating the energy per site 1 < n < N as 

E n = t^pI + ~V(q n+ i - q n ) + ^V(q n - q n - X ) (15) 

where V(x) — ^x 2 + ^x 4 , thereby defining the function 

C Q {t) = N ^jr 1 i . (16) 

Since, Co is of order one if Ei — E/N at each site i of the chain and of order N 
if the energy is localized at only one site, this function can serve as an efficient 
indicator of energy localization in the chain. 

In their experiments, Cretegny et al. used TV = 128, (3 — 0.1, E > E u , 
where E u is the lowest destabilization energy of the 7r-mode and plotted Co 
versus t. Distributing evenly the energy among all sites at t = according to 
the 7r-mode pattern, they observed that Co initially grows to relatively high 
values, indicating that the energy localizes at a few sites. After a certain time, 
however, Co reaches a maximum and decreases towards an analytically derived 



asymptotic value Co — 1.795 [35[, which is associated with the breakdown of 
the chaotic breather and the onset of energy equipartition in the chain. 

We have performed the same study of the 7r-mode at a lower energy, E = 
0.768 > E u w 0.0257, (/3 = 1 and N = 128 in the present work) and have 
been able to relate the g-Gaussian distributions of our QSS to the lifetime of 
chaotic breathers reported in [35| as follows: In panel (a) of Fig. [31 we plot 
Co as a function of time and verify indeed that it grows over a long interval 
(t ~ 1.8 x 10 8 ) after which it starts to decrease and eventually tends to the 
asymptotic value Co — 1.795 associated with energy equipartition. In panel (b) 
of the same figure we plot the evolution of the four biggest Lyapunov exponents 
and observe that initially they decrease towards zero until about t ~ 4x 10 7 when 
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Figure 3: Panel (a) is a plot of Co as a function of time for the unstable (E = 0.768 > 
E u fa 0.0257) 7r-mode with /3 = 1 and N = 128. Panel (b) is a log-log plot of the four 
biggest Lyapunov exponents as a function of time for the same parameters as in panel (a). 
Repeating the calculation of the exponents for time windows of different size preserves the 
qualitative nature of the plot and changes only the smaller values of the exponents by one 
order of magnitude. 



they start to increase towards positive values indicating the unstable character 
of the 7r-mode. 

We now present in Fig. S]the instantaneous energies E n of all N — 128 sites 
at two different times. In the first panel, at time t = 10 7 , where Co(i) attains 
a maximum (see Fig. Ufa)), the energy is clearly localized at very few sites, 
demonstrating the occurrence of a chaotic breather. At later times, however, 
e.g. t — 6x 10 8 , panel (b) shows that this breather has collapsed, as the system 
tends to equipartition. This is consistent with the fact that the pdf describing 
the statistics at t = 10 7 is well-approximated by a g-Gaussian, with q w 2.6, 
shown in panel (c), while at t = 10 s , as we know from Fig. QJc), this pdf is 
already very close to a true Gaussian. 

Finally, in Fig. EJd) we present estimates of the q values one obtains, when 
computing chaotic orbits near the 7r-mode at this energy density e = E/N — 
0.006. Remarkably enough, even though these values have an error bar of about 
±10 percent due to the different statistical parameters (M,N^ C ) used in the 
computation, they exhibit a clear tendency to fall closer to 1 for t > 10 7 , where 
energy equipartition is expected to occur. 

It is indeed a hard and open problem to determine exactly how equipartition 
times Teq scale with the energy density e = E/N and other parameters (like /3), 
particularly in the thermodynamic limit, even in 1-dimensional FPU lattices. 
Although it is a question that has long been studied in the literature (35l . [Ici - tEH , 
the precise scaling exponents by which T e q depends on e, j3, etc. are not yet 
precisely known and the controversy is continuing to this very day. 

It would be very interesting if our q values could help in this direction. In 
fact, carrying out more careful calculations, as in Fig. HJd), at other values of 
the specific energy, e.g. e = 0.04 (see Fig. [5]) and e = 0.2, we found q plots 
that exhibited a clear decrease to values close to 1, at T e q ~ 7.5 x 10 5 and 
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Figure 4: In panel (a) at t = 10 7 , near the maximum of Co(t) (see Fig. Efa)), we see a chaotic 
breather. In panel (b) at t = 6 X 10 s , this breather has collapsed and the system has reached 
a state whose distribution is very close to a Gaussian, as implied already by the pdf shown in 
Fig. fTT c) at t = 10 s . Note the scale difference in the vertical axes. Panel (c) at t = 10 7 shows 
that the sum distribution, near the maximum of the chaotic breather, is still quite close to 
a tj-Gaussian with q 2.6. Panel (d) presents an estimate of the q index at different times, 
which shows that its values on the average fall significantly closer to 1 for t > 10 7 . 



T e q pa 7.5 x 10 4 respectively, approximately where the corresponding chaotic 
breathers collapse. Still, even though our results are consistent with what is 
known in the literature, the limited accuracy of our approach does not allow 
us to say something meaningful about scaling laws, as the system tends to 
equilibrium. In any case, the quantitative usefulness of our q estimates is a very 
interesting topic, which we are currently exploring and intend to describe in 
more details in a future publication. 

In conclusion, comparing Figs. [TJ |3][a) and HI we deduce that while the 
chaotic breather still exists, a QSS is observed fitted by q-Gaussians with q well 
above unity, as seen in Fig. [4jc) . However, as the chaotic breather breaks down 
for t > 10 s and energy equipartition is reached (see Fig. Sfb)), q 1 and 
g-Gaussian distributions rapidly converge to Gaussians (see Fig. [Tfc)) in full 
agreement with what is expected from BG Statistical Mechanics. 
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Figure 5: Plot of q values as a function of time for energy density e = 0.04, i.e E = 5.12 (green 
curve) for the 7r-mode of an FPU periodic chain with TV = 128 particles and j3 = 1. The 
blue and red curves show error bars in the form of plus and minus one standard deviation. In 
agreement with other studies the transition to values close to q = 1 occurs at Teq ~ 7.5 X 10 5 . 



5. FPU SPOl and SP02 modes under fixed boundary conditions 

Let us now pass to the next part of our study and examine the chaotic 
dynamics near NNMs of the FPU system under fixed boundary conditions. In 
particular, we first study pdfs of sums of chaotic orbit components near a NNM 
we call the SPOl mode (see 0, [H^]), which keeps one particle fixed for every 
two adjacent particles oscillating with opposite phase. 

More specifically, we consider again the FPU-/3 1-dimensional lattice de- 
scribed by the Hamiltonian (|10|) and impose fixed boundary conditions 

q (t) =q N +l(t) =0,Vt (17) 

with N odd. The SPOl mode is defined by 

N — 1 

q 2j (t)=0, q2j-i(t) = -q 2 j+i(t), j = 1,..., — — . (18) 

As shown in Fig. [6l the chaotic region very close to this solution (when it 
has just become unstable) appears for a long time isolated in phase space from 
other chaotic domains. In fact, one finds several such domains, embedded one 
into each other. For example, in the case N — 5 with j3 = 1.04 a "figure eight" 
chaotic region appears in blue traced by an orbit starting at a distance about 
1.192 x 10~ 7 from the SPOl mode. This is easily seen on the surface of section 
(lii Pi) °f Fig- HI computed at times when q 3 = (and p$ > 0) at the energy 
E = 7.4. Even though the SPOl mode is unstable (depicted as the saddle point 
for q 3 — 0) orbits starting sufficiently nearby remain in its vicinity for very long 
times, forming eventually the thin "figure eight" at the center of the figure. 

Starting, however, at points a little further away (e.g. at an initial condition 
located at a distance about 1.086 x 10~ 2 from the SPOl mode), a more extended 
chaotic region is observed, plotted by green points, which still resembles a "figure 
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eight" . Choosing even more distant initial conditions (e.g. one that is located 
at a distance about 3.421 x 10 _1 from the SPOl mode), a large scale chaotic 
region plotted by red points becomes evident in Fig. |B1 

It is, therefore, reasonable to regard these different dynamical behaviours 
near the SPOl mode as QSS and characterize them by pdfs of sum distributions, 
as explained in Section [2] The idea behind this is that orbits starting initially 
at the immediate vicinity of the unstable SPOl mode behave differently than 
those lying further away, since the latter orbits have the ability to explore more 
uniformly ergodic parts of the energy manifold. We thus expect to find, for the 
nearby regions of SPOl, q-Gaussian-like distributions with 1 < q < 3, while 
for orbits starting sufficiently far the distributions are expected to be Gaussians 
with q — > 1. 

2.8 
2.6 
2.4 

a 22 
2 

1.8 

1.6 

-0.2 -0.1 0.1 0.2 

11 

Figure 6: The "figure eight" chaotic region (blue colour) for an initial condition at a distance 
about 1.192 X 10 — 7 from the unstable SPOl mode (depicted as the saddle point when 53 = 
and P3 > 0), a slightly more extended "figure eight" region (green points) for an initial 
condition a little further away (ri 1.086 X 10 -2 ) and a large scale chaotic region (red points) 
for an initial condition even more distant (fs 3.421 X 10 _1 ) on the surface of section (<ji,pi) 
computed at times when 93 = 0. We have integrated our orbits up to tf = 10 5 using E = 7.4, 
N = 5 and /3 = 1.04. 

To test the validity of these ideas, we have chosen the quantity 

V(t) = qi(t) + q 3 (t) (19) 

as our observable, which is exactly equal to zero at the SPOl orbit (see Eq. (|18p). 
In fact, rj(t) remains very close to zero in the energy interval where the SPOl 
mode is stable and becomes non-zero at energies beyond the first destabilization 
energy of the mode. Following what we have presented in Section [21 we now 
study the three different initial conditions located in the neighborhood of the 
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unstable SP01 mode, seen in Fig. [51 In particular, in panel (a) of Fig. we 
see the surface of section created by the trajectory starting at a distance about 
1.192 x 10~ 7 from it and integrated up to tj? = 10 5 while in the following two 
panels we see the same surface of section computed for final integration times 
of tc = 10 7 and — 10 8 respectively. The parameters are the same as in Fig. [51 
It is obvious that, as the integration time increases, even orbits which start 
close to the unstable SPOl mode, eventually wander over a more extended 
part of phase space, covering gradually all of the energy manifold when the 
integration time is sufficiently large (e.g. if = 10 8 ). This can also be explained 
by the Lyapunov exponents seen in panel (d). While for a fairly long time 
interval (up to approximately t — 10 7 ) they decrease towards zero, they suddenly 
jump to much higher values after t « 10 7 , indicating that the orbit has entered 
into a wider chaotic domain of phase space. 
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Figure 7: Panel (a) shows the (<?i,pi) surface of section of an orbit integrated up to t£ = 10 5 
and starting at a distance about 1.192 X 10 — 7 from the unstable SPOl for N = 5 and /3 = 1.04 
at the energy E = 7.4. Panel (b) is the same as panel (a) but for = 10 7 . Panel (c) is the same 
as in panel (a) but for if = 10 s . Panel (d) shows the corresponding four biggest Lyapunov 
exponents. Plots in linear- log scale of numerical (red curve), of g-Gaussian (green curve) 
and Gaussian (blue curve) distributions for the initial condition corresponding to panels (a) 
to (d). In panel (e) we use 4r = 10 5 , iVj c = 10 4 and M = 10. In this case the numerical 
fitting with a g-Gaussian gives q Rj 2.785 with x 2 ~ 0.000 31. In panel (f) integration time 
is t£ = 10 7 , 7Vj c = 10 5 and use M = 1000 terms. Here, the numerical fitting gives q fa 2.483 
with x 2 ~ 0.000 47. In panel (g) we set t { = 10 s , N ic = 10 5 and M = 1000 and obtain a 
distribution much closer to a Gaussian (q Ri 1.05), corresponding to the chaotic orbit shown 
in panel (c). 
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An important question arises now: Are these different behaviours reflected 
in the statistical distributions associated with these trajectories computed for 
successively longer times? The answer to this question is presented in the three 
lower panels of Fig. [7] In particular, in panel (e) we computed in linear-log 
scale the numerical (red curve), q-Gaussian (green curve) and Gaussian (blue 
curve) distributions for the initial condition located closest to the SPOl, using 
if = 10 5 , iVj c = 10 4 time windows and M = 10 terms in the computations of 
the sums. In this case, a fitting with the q-Gaussian ([7]) gives q s» 2.785 with 
X 2 ~ 0.000 31. This distribution corresponds to the surface of section shown in 
panel (a). If we now increase if by two orders of magnitude (see panel (f)) using 
iVj c = 10 5 , M = 1000 and perform the same kind of fit we get q ps 2.483 with 
X 2 0.000 47. 

It is important to emphasize that the lower parts (tails) of the red distribu- 
tion of panel (f) are not fitted well by the q-Gaussian. This suggests that by 
increasing the integration time, the initial g-Gaussian takes a transient form, 
which may very well be approximated by other types of functions like Eq. (|14[) . 
This distribution corresponds to the surface of section of panel (b) . By increas- 
ing the final integration time further to if = 10 8 and using N^ c = 10 5 and 
M = 1000 terms, we observe that the red curve of panel (g) is very close to a 
Gaussian (q sa 1.05), characterizing the chaotic regime plotted in the surface of 
section of panel (c). 

Let us now perform the same kind of study for the second initial condition 
of Fig. |6l starting at a distance about 1.086 x 10~ 2 from the SPOl periodic 
orbit. In this case we also obtain very similar results as in Fig. [7l Again, the 
lower parts (tails) of the distribution deviate clearly from a q-Gaussian shape. 
In fact, for if = 10 8 , they are even closer to a Gaussian, as the corresponding 
orbit covers a much larger part of the available energy manifold. 




Figure 8: Panel (a) is the (gi,pi) surface of section of an orbit in an initial distance about 
3.421 X 10 _1 from the unstable SPOl orbit integrated up to = 10 5 . In panel (b) a plot of 
the data in linear-log scale shows that it is very well-fitted by a Gaussian (blue curve) with 
q es 1.0. The final integration time is = 10 5 using 2Vj c = 10 5 time windows and M = 1000 
terms in the computation of the sums. 
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Figure 9: Panel (a) is the (q±,pi) surface of section of an orbit integrated up to t£ = 10 10 and 
starting initially at a distance about 1.418 X 10 -3 from the unstable SP02 mode. Panel (b) 
shows the corresponding four biggest Lyapunov exponents. Panel (c) presents in linear-log 
scale the numerical (red curve), g-Gaussian (green curve) and Gaussian (blue curve) distri- 
butions. The final integration time is tf = 10 10 using N\~ = 333 333 333 time windows and 
M = 10 terms in the sums. Here, the numerical fitting with a q-Gaussian gives q r& 2.769 with 
X 2 ~ 4.44 X 10 — 5 , but panel (d) shows that the red distribution of panel (c) is better fitted 
by the P function of Eq. (X3J for ai « 0.006, a q fa 170 and q fa 2.82 with \ 2 fa 2.06 X 10 -6 
(green curve). 



Finally, we have carried out the same analysis for the third initial condition 
of Fig. [6] located at a distance about 3.421 x 10 _1 from the SPOl mode. This 
corresponds to the orbit yielding the chaotic region shown in Fig. HJa). Note 
that despite the large size of this region compared to the ones shown in Fig. 
[Tfa) and (b) it is deeply embedded in the chaotic domain of Fig. Etc). Due to the 
presence of "strong chaos" here, we find that even for small integration times 
(tf = 10 5 ), the distribution is very well approximated by the Gaussian shown in 
Fig.Efb). 

In order to study further the persistence of QSS as time increases, we turned 
to another nonlinear mode of the FPU Hamiltonian with fixed boundary con- 



ditions called the SP02 mode [53j. This is a NNM which keeps every third 



particle fixed, while the two in between move in exact out of phase motion. 
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What is interesting about this NNM is that it becomes unstable at much lower 
energies (i.e. E^/N cx 7V~ 2 ) compared to SPOl {E\/N oc iV" 1 ) Q, much like 
the low k = 1,2,3^. . mode periodic orbits connected with the breakdown of 
FPU recurrences [UIH^. Thus, we expect that near SP02, orbits will be more 
weakly chaotic than SPOl and hence QSS are expected to persist for longer 
times. This is exactly what happens. As Fig. [S] clearly shows, the dynamics in 
a close vicinity of SP02 has the features of what we might call "edge of chaos" : 
Orbits wander in a regime of very small (positive) Lyapunov exponents, tracing 
a kind of "banana" shaped region much different than the "figure eight" we had 
observed near SPOl, while their pdfs (up to if = 10 10 ), converge to a function 
that is close to a g-Gaussian, never deviating towards a Gaussian, as we have 
seen in the QSS of other FPU systems studied in this paper. 

More specifically, let us take N = 5, (3 = 1, E — 0.5 and choose an orbit 
located initially at a distance about 1.418 x 10 -3 from the SP02 solution, which 
has just turned unstable (at E\ ss 0.4776). As we can see in Fig. HJa), the 
dynamics here yields banana- like orbits at least up to if = 10 10 . The weakly 
chaotic nature of the motion in this domain is plainly depicted in panel (b) of 
this figure, where we have plotted the four positive Lyapunov exponents up to 
if = 10 10 . Note that, although they all decrease towards zero for a very long 
time, at about if > 10 9 , the largest one of them tends to converge to a very 
small value (about 10~ 8 ), indicating that the orbit is chaotically sticking to an 
"edge of chaos" region around SP02. 

In panel (c) of this figure, we have plotted the corresponding pdf at time 
if = 10 ln (the function does not change after if = 10 7 ). What we discover is 
an extremely long-lasting QSS, whose distribution is well-fitted by a g-Gaussian 
with q pa 2.769 and \ 2 ~ 4.44 x 10 -5 . The "legs" of the distribution away 
from the center deviate from the g-Gaussian shape, but remain very far from 
the Gaussian function also plotted in the figure by blue colour. We have also 
performed a similar fitting of our data with the function (|14[) . as in the case 
of panel (d) of Fig. [T] What we find here is that the numerical distribution 
(red curve) of panel (c) is more accurately approximated by Eq. (|14[) where 
ai fa 0.006, a q pa 170 and q pa 2.82 with X 2 ~ 2.06 x 10~ 6 , in contrast to the 
X 2 pa 4.44 x 10~ 5 obtained by fitting the same distribution by a g-Gaussian with 
q pa 2.769 (see Fig. Efa)). 

Thus, we conclude that in sufficiently thin chaotic layers of multi-dimensional 
Hamiltonian systems with very small positive Lyapunov exponents it is possible 
to find non-Gaussian QSS that persist for quite long times as in the SP02 
case. Our numerical evidence suggests that in these regimes chaotic orbits stick 
for long time intervals to a complex network of islands, where their statistics 
is well approximated by distributions of the g-Gaussian type connected with 
Nonextensive Statistical Mechanics. 

It is now interesting to study the dynamics near the unstable SPOl and 
SP02 modes using an analysis similar to the one carried out for the 7r-mode in 
Section [4l following [35j |. In particular, focusing on small perturbations of both 
modes under fixed boundary conditions, we wish to see how Co (see Eq. (fl"6|) ) 
displays the evolution towards energy equipartition and whether this transition 
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will again be preceded by the appearance of chaotic breathers. 

To find out, we used N = 129 particles with /? = 1.04 and E = 1.5 for the 
SPOl mode and considered a neighboring orbit starting at an initial distance 
2.22 x 1CT 7 from it. For these parameters, ^ P01 ~ 1.05226 [Hj]. As we see 
in Fig. [TU1 Co grows on the average (without a clear maximum) and settles 
down to a value near 1.75 indicating that the system reaches equipartition at 
about t « 10 6 , when the four biggest Lyapunov exponents begin to converge to 
their final values (see panel (b) of Fig. [TU]) . Interestingly, during this period, 
the corresponding distribution in panel (c) (using the observable rj — (?(64) + 
g(65) + <?(66)) is well fitted by a q- Gaussian with q 2.843 and \ 2 ~ 0.0003 at 
t = 4x 10 5 . For longer times, distributions quickly tend to Gaussians as shown 
in panel (d) for t = 2 x 10 6 . Thus, even though (unlike the 7r-mode) no chaotic 
breather is observed here, during the time Co takes to relax to its limiting value, 
sum distributions are well approximated by q-Gaussians and only converge to 
Gaussians when energy equipartition has occurred. 

Let us repeat this process now for an orbit starting at an initial distance of 
4.09 x 10~ 7 from the SP02 mode, with N = 128, = 1 and E = 0.1, above 
the first destabilization threshold of SP02 (i.e. #SP° 2 « 0.01279) j^. As 
we observe in Fig. 111! Co also does not exhibit a distinct high maximum here, 
but grows on the average more slowly than in the SPOl case. Thus, it takes 
longer to approach its limit, indicating that the system reaches equipartition 
at t s» 4 x 10 8 , where the four biggest Lyapunov exponents cease to decrease 
towards zero and tend to small positive values (see panel (b) of Fig. [TTj) . The 
corresponding sum distribution shown in panel (c) at t = 10 6 is still well fitted 
by a q-Gaussian with a high q = 1.943 value, while we need to increase the time 
to 5 x 10 8 to see pdfs that are much closer to a Gaussian (q w 1.0), as in panel 
(d), indicating that this is a lower bound for energy equipartition. 



6. q-Gaussian distributions for a small microplasma system 

Let us turn our attention, finally, to a Hamiltonian system, which is very 
different than the FPU systems of the previous sections. In particular, we shall 
work with a system of few degrees of freedom characterized by long range in- 
teractions of the Coulomb type. This will allow us to study statistically the 
dynamics of its transition from a crystalline-like to a liquid-like phase (the so- 
called "melting transition") [37|,|38| at small energies, as well as identify a tran- 
sition from the liquid-like to a gas-like phase, as the total energy is constantly 
increased [HB|. Before presenting our results, we first describe the model and its 
various geometries. 

Our microplasma system consists of N ions of equal mass m = 1 and electric 
charge Q moving in a Penning trap with electrostatic potential 37, 55| 



2z^ — — v 

$(x,y,z) = V / (20) 
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Figure 10: Panel (a) is a plot of Co (red curve) as a function of time for a perturbation of the 
unstable (E = 1.5 > £;? P01 ~ 1.05226) SPOl mode with /S = 1.04 and N = 129 at initial 
distance 2.22 X 10~ 7 . The green curve corresponds to the time average of the red curve. Panel 
(b) is a log-log scale plot of the four biggest Lyapunov exponents as a function of time. Plots in 
linear-log scale of numerical (red curve), g-Gaussian (green curve) and Gaussian (blue curve) 
distributions: In panel (c) for final integration time ([ = 4 X 10 5 , N^ c = 10 4 time windows 
and M = 20 terms in the sums, the numerical fitting with a q-Gaussian gives q fa 2.843 with 
X 2 fa 0.0003, while panel (d) corresponds to t f = 2 X 10 6 , N ic = 5 X 10 4 and M = 20. Here, 
the numerical fitting gives q 1.564 with \ 2 ~ 0.00014. 



and a constant magnetic field in the z direction, whose vector potential is 

1 



A(x,y,z) = ~(-By,Bx,0). 
The system is described by the Hamiltonian 

« = E(^lw-«A(r i )] 2 + O*(r0}+ £ 

i=l I J l<i<j<N 



Q 2 



-iire rij 



(21) 



(22) 



where i"i is the position of the zth ion, is the Euclidean distance between zth 
and j'th ions and £q is the vacuum permittivity. In the Penning trap, the ions 
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Figure 11: Panel (a) shows a plot of Co (red curve) as a function of time for a perturbation of 
the unstable (E = 0.1 > E^ P ° 2 ks 0.01279) SP02 mode with = 1 and TV = 128 at initial 
distance 4.09 X 10" 7 . The green curve corresponds to the time average of the red curve. Panel 
(b) presents in log-log scale the four biggest Lyapunov exponents as a function of time for the 
same parameters and initial condition as in panel (a). In panel (c) we plot in linear- log scale 
of numerical (red curve), q-Gaussian (green curve) and Gaussian (blue curve) distributions 
for the same initial condition and parameters as in panel (a) for if = 10 6 , iVj = 2.5 X 10 4 
and M = 20. Here, the numerical fitting gives q ss 1.943 with x 2 ~ 0.00035, which is further 
from a Gaussian than the corresponding pdf for SPOl shown in Fig. llOf d). Panel (d) shows 
that the pdfs are much closer to Gaussian (q rj 1.0) at if = 10 8 , when the system has nearly 
reached the equipartition state. 



are subjected to a harmonic confinement in the z direction with frequency 



4QV 



m(rg + 2z 2 ) 



(23) 



while, in the perpendicular direction (due to the cyclotron motion) , they rotate 
with frequency w c = QB/m. Thus, in a frame rotating about the z axis with 
Larmor frequency u-^ — w c /2, the ions are subjected to an overall harmonic 
potential with frequency uj x = u y = (wc/4 — wl/2) 1 / 2 in the direction perpen- 
dicular to the magnetic field. In the rescaled time r = cj c i, position R = r/a 
and energy H = r H/{muiQa 2 ) with a = [Q 2 /(47reomwc)] 1 ^ 3 , the microplasma 
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Hamiltonian (|2"21 takes the form 



N N 



» -5 £ p = + ml - t) (X - + + T z ?] + E i = £ < 24 » 

i—1 i—1 i<Cj J 

where E is the total constant energy, = (Xi, Y^Zi) and Pi = [Px t i Pyi i ) 
are the positions and momenta in R 3 respectively of the N ions, Rij is the 
Euclidean distance between different ions i,j given by 

Rij = [(X, - Xjf + (Yi - Yjf + (Zi - Zjf\^ (25) 

and finally 7 = lu z /loc- 

The ions perform bounded motion under the condition that 

< l7l < ^= (26) 

and the Penning trap is called prolate if < [7] < l/v6, isotropic if I7I = l/v6 
and oblate if 1/ V6 < I7I < 1/ Thus, the motion is quasi 1-dimensional in the 
limit 7 — > and quasi 2-dimensional in the limit 7 — > 1/V2. The Z direction 
is a symmetry axis and hence, the Z component of the angular momentum 

= J2iLi XiPyt — YiPxi is conserved, being a second integral of the motion. 
We set, from here on, the angular momentum equal to zero (i.e. = 0) and 
study the motion in the Larmor rotating frame. 

In [37I ] , the authors demonstrate the occurrence of a dynamical regime change 
in a system (|24|) composed of N = 5 ions and confined in a prolate quasi 1- 
dimcnsional configuration of 7 = 0.07. More specifically, in the lower energy 
regime, a transition from crystalline-like to liquid-like behaviour is observed, 
called the "melting phase". It was first shown that this melting of the crys- 
tal is not associated with a sharp increase of the temperature at some critical 
energy, as might have been expected in first sight. Furthermore, the positive 
Lyapunov exponents are maximal at energies much higher than the "melting" 
regime. Thus, it appears that no clear macroscopic methodology is available for 
identifying and studying this melting process in detail. 

For these reasons, the Smaller Alignment Index (SALI) method (56l - [58l was 
used by the authors to study the local microscopic dynamics in detail [37j. It 
was discovered that there exists an energy range of weakly chaotic behaviour, 
i.e. E G AE 1Ti ^ = (2,2.5), where the positive Lyapunov exponents are very 
small and the SALI exhibits a stair-like decay to zero with varying decay rates. 
This suggests the presence of long-lived "sticky" orbits executing a multi-stage 
diffusion process near the boundaries of resonance islands 15]] . Thus, it was 
concluded that it is in this energy interval that the melting transition of the 
system (f2"4"| occurs. 

Motivated by the above results and pursuing this investigation further, we 
decided to restrict ourselves to the same 7 and N, at which the minimum en- 
ergy for ions to start moving appreciably about their equilibria positions is 
Eq R3 1.8922. Our first aim was to study the melting transition of the system 
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as the energy of the Hamiltonian (pM]l increases above E , using pdfs associ- 
ated with chaotic trajectories to relate microscopic to macroscopic observables. 
Based on the results of the previous sections, we expected to find g-Gaussian 
approximations with 1 < q < 3 in the vicinity of AE m i , where the positive Lya- 
punov exponents A.; , i = 1, . . . , 3iV are quite small compared to their maximum 
value of Ai 



37] 



0.0558 attained atfiw 5.95 
In Fig. [121 we show the results of our study for the interval (So, 10]. We 
chose as an observable the quantity rj(t) — X\(t), i.e. the first component of 
the position Ri of the first ion and set A/j c = 2 x 10 4 , M = 1000 and total 
integration time ij = 2x 10 7 . Just as expected, in the energy range AS m ^ of 
the melting transition, the values of the entropic parameter of the associated 
g-Gaussian pdf were well above q = 1, indicating that the statistics is certainly 
not Gaussian. In fact, the detected g-Gaussian shape persists up to E w 4, 
implying perhaps that the energy interval of the melting transition is slightly 
bigger than originally estimated in [37 1. 
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Figure 12: Plot in log-linear scale of the q-entropic parameter (red colour) as a function of 
the energy E of the microplasma J24D for 7 = 0.07 (prolate trap) and N = 5. In this plot, we 
have used N- lc = 2 X 10 4 , M = 1000 and have integrated every orbit up to = 2 X 10 7 . We 
have also plotted the line at q = 1 (denoted by green) for comparison with the Gaussian case. 



Next, we proceeded to study the second dynamical state transition of system 
(|24]l. from a liquid- like to a gas-like regime, where the system evolves in a 
"weakly chaotic" way over an energy range where the largest positive Lyapunov 
exponents decrease towards zero following Eq. ([27| , as shown in Fig. ITBT a) . We 
thus move to higher energies and concentrate on the passage of our system from 
a strongly chaotic regime at energies where the Lyapunov exponents attain their 
maximum (i.e. for E > 6) to energies where the motion is much less chaotic. 
In that regime, the highly chaotic but strongly correlated state of frequent 
inter-particle close encounters ("collisions") is replaced by considerably fewer 
inter-particle energy/momenta exchanges and consequently much weaker chaos. 

In 55j it is demonstrated that if all TV 2 inter-particle interaction terms of the 



Coulomb part of Hamiltonian (l24l) at distances Rij contribute to the maximal 
Lyapunov exponent Ai, one derives, for sufficiently large values of T — > 00 (and 
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therefore of E — > oo), that 



(27) 



where T is related to the temperature of the system. This formula explains the 
asymptotic power law decay of the biggest Lyapunov exponent observed in panel 
(a) of Fig. Q21 It can also be seen that the second largest Lyapunov exponent 
A2 (plotted with green colour) obeys a similar formula as a function of E. 

It is, therefore, very interesting that the q indices of the corresponding pdfs 
remain well above unity for E £ (30, 200), as shown in Fig. [T3](b) . This suggests 
that in this range there is a significant dynamical change of the system to 
a weaker form of chaos, associated with small positive Lyapunov exponents. 
Most probably, the N 2 inter-particle terms do not contribute significantly to 
the transition between the fully developed chaotic regime of correlated motion 
and the more ordered dynamics of a gas-like system. Thus, we conclude that the 
energy increase drives our system to a more and more organized state, favoring 
the kinetic over the Coulomb part, where very few close encounters between the 
particles occur, just like the case of an "ideal gas" . 
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Figure 13: Panel (a) shows the two biggest positive Lyapunov exponents (Ai with red and 
A2 with green colour) as a function of the energy E in log-log scales. We observe that, after 
attaining their peak values at E 5.8, both of them decay to zero according to formula H27II . 
Panel (b) is a plot in log-linear scales of the q-entropic parameter (red colour) as a function 
of the energy E of the microplasma Hamiltonian I I24H for 7 = 0.07 and TV = 5. We have used 
N ic = 2 X 10 4 , M = 1000, i[ = 2x 10 7 and also plotted the line at q = 1 (denoted by green) 
for reference to the entropic parameter of the Gaussian distribution. 



7. Conclusions 

In this paper, we have numerically constructed pdfs of rescaled sums of 
observables derived from chaotic trajectories of multi-dimensional Hamiltonian 
systems. Assuming that these observables behave as independent random vari- 
ables, we sought to determine their statistics in phase space regions of "weak 
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chaos" in the context of the Central Limit Theorem. More specifically, we fo- 
cused on the neighborhood of unstable periodic orbits of various examples of 
FPU-/? particle chains and a microplasma system undergoing phase changes 
from crystalline to liquid and from liquid to gas-like dynamics, as its total en- 
ergy increases. 

For the FPU systems, we discovered that when chaotic orbits are restricted 
within "small size" domains near NNM which have just turned unstable, their 
pdfs are well approximated by g-Gaussian functions (1 < q < 3) for long times 
(typically up to t ps 10 6 ). These pdfs frequently represent QSS and tend to 
Gaussians (q = 1) for longer integration times, as the orbits diffuse into do- 
mains of "strong chaos" characterized by large (positive) Lyapunov exponents. 
There are also cases, however, where the initial conditions are so close to an 
"edge of chaos" regime that the pdfs converge to a smooth function, which is 
well- approximated by different analytical expressions, valid in crossover regimes 
between g-Gaussians and Gaussians, as discovered recently also by other re- 
searchers in very different systems. 

Thus, it is important to point out that g-Gaussians offer a family of func- 
tions, parametrized by a single index q, that can successfully approximate sum 
distributions of QSS evolving in complicated networks of chaotic motion of 
multi-dimensional Hamiltonian systems. Indeed, g-Gaussians may represent the 
"best" choice for uniformly fitting statistical data in cases where the passage 
from an "edge of chaos" regime to a domain of uniformly spread chaos is identi- 
fied by a sequence of q values starting well above 1 and tending to unity as time 
increases. With regard to the relation of g-Gaussians to certain important phys- 
ical properties of FPU systems, such as energy equipartition, we have observed 
that these functions accurately represent long-lived transient states, whose q is 
significantly larger than 1 and tend to Gaussians with q falling quickly to unity, 
as the energy begins to be equally shared by all degrees of freedom. 

In the case of the microplasma Hamiltonian, the approximation of our dis- 
tributions by g-Gaussians allowed us to identify two distinct energy regimes of 
"weakly" chaotic behaviour (and small positive Lyapunov exponents): A low 
energy interval, E € [2, 4], over which the system undergoes a melting transition 
and a much longer range, E £ [30, 200], where it passes from liquid-like to gas- 
like dynamics. This is accomplished by observing over which energy domains 
the g-index of the distributions attains values significantly higher than unity, 
before returning to the q — 1 value of Gaussians characterising "strong chaos" . 

Other authors have also studied QSS in conservative systems from the point 
of view of Nonextensive Statistical Mechanics 2l|, |22j , but do not compare sum 
distributions of their orbits, as we have done here. Rather, they calculate a 
quantity T(t) called "dynamical temperature", which measures the total "an- 
gular momentum" of their system of N coupled standard maps. They show 
that: (i) in the N = 1 case, T(t) grows smoothly to a ^Qgg(a) value, which 
depends on the nonlinearity parameter a and is smaller, for a = O(l), than 
the TgQ = 1/12 expected for a BG system, and (ii) for N > 2-dimensional 
maps, T(t) — > TgQ typically over 10 4 iterations for 0.3 < a < 0.5. Their initial 
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conditions are chosen randomly within a thin layer that includes chaotic as well 
as regular orbits, thus it is not clear how to interpret their QSS and verify their 
results in more complicated systems. For example, the ?Qgg(a) < Tqq values 
for a single standard map may correspond to cases where most orbits remain 
within islands of regular motion. 

On the other hand, in our study of FPU Hamiltonians, we have discovered 
QSS, which either diffuse to larger domains passing from a g-Gaussian to a Gaus- 
sian state, or may converge to a non-Gaussian type of pdf when the Lyapunov 
exponents are vanishingly small. It would be interesting, therefore, in a future 
publication, to define an analogous "temperature" quantity for our Hamiltonian 
systems, as was done in 2l], 22 1 and study the times it needs to reach the BG 
state, as a function of the width of the layer of initial conditions. 

Interestingly, certain very recent results on low-dimensional area-preserving 
maps [36| are in very good agreement with what we have found here for multi- 
dimensional Hamiltonian systems. These results suggest that the presence of 
families of islands of stability in phase space and the associated stickiness phe- 
nomena in weakly chaotic regions of conservative systems are responsible for 
QSS approximated by g-Gaussian distributions. However, as soon as dissipa- 
tive effects are included and strange attractors appear, chaotic motion becomes 
more uniform and sum distributions converge rapidly to Gaussians. 

We believe, therefore, that the realm of conservative systems, represented 
by Hamiltonian flows or symplectic maps, is well suited for discovering com- 
plex metastable phenomena occurring in "weakly chaotic" domains. Even if 
the system is multi-dimensional, the dynamics near the boundaries of resonance 
islands and small positive values of Lyapunov exponents near these boundaries 
frequently yield long-lived states with nonextensive statistics. Thus, we suggest 
that it may be possible to apply these ideas to some very slow diffusive phenom- 
ena recently observed in 1-dimensional disordered chains T^- 18| , where it is not 



yet known if, in the limit of t — > oo, diffusion will extend to particles further 
and further away or will be hindered by the boundary of some high-dimensional 
KAM torus, which is believed to exist 15911 - 
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